function C = contingency_matrix (c1, c2, rad)

	if (isscalar(rad))
		rad = [rad rad];
	end

	C = zeros(c1.topol(1)*c1.topol(2), c2.topol(1)*c2.topol(2));
	for i = 1:c1.topol(1)*c1.topol(2)
		for j = 1:c2.topol(1)*c2.topol(2)
			data1 = get_data_from_neuron(c1.bmus, i, rad(1), c1.topol);
			data2 = get_data_from_neuron(c2.bmus, j, rad(2), c2.topol);
			C(i, j) = length(intersect(data1, data2));
		end
	end